In Metropolitan-Council/planting.shade: Tree planting and preservation prioritization tool

block group issues

the current watermask in gee is only for metcouncil area; this dataset probably works for rest of state: https://gisdata.mn.gov/dataset/water-national-hydrography-data (but need something from WI in there too)

would want to focus only on residential area or city boundaries; ag land has different purposes/uses/challenges/sustainability goals.

very slow when showing block groups across entire state

-- options: focus on statistical areas; some of this won't be super useful either since ag land focus on areas where population density is above some threshold

knitr::opts_chunk$set(echo = FALSE, message = FALSE, warning = FALSE,
                      cache = TRUE,
                      cache.path = "cache/")
library(tidyverse)
library(tigris)
library(sf)
library(cowplot)
library(leaflet)
st_erase = function(x, y) st_difference(x, st_union(st_combine(y)))
`%not_in%` <- Negate(`%in%`)

df1 <- bg_growingshade_main %>%
  filter(variable %in% c("pbipoc", "canopy_percent", "mdhhincnow", "avg_temp", "ndvi", "pop_density")) %>%
        select(tract_string, variable, raw_value) %>%
        pivot_wider(names_from = variable, values_from = raw_value)

df <- df1 %>%
        select(canopy_percent, mdhhincnow, pbipoc) %>%
        pivot_longer(names_to = "names", values_to = "raw_value", -c(canopy_percent)) %>%
        mutate(raw_value = if_else(names == "pbipoc", raw_value * 100, raw_value))

fig_equity <-
        ggplot(aes(x = raw_value, y = canopy_percent), data = df) +
        geom_point(col = "grey40", alpha = .2, data = filter(df), na.rm = TRUE, size = .7) +
        geom_smooth( # method = "lm",
          method = "gam", formula = y ~ s(x, bs = "cs"),
          fill = NA, col = councilR::colors$councilBlue, na.rm = TRUE
        ) +
        councilR::council_theme()  +
        facet_wrap(~names,
          scales = "free", nrow = 1, strip.position = "bottom",
          labeller = as_labeller(c(pbipoc = "Population identifying as\nperson of color (%)", mdhhincnow = "Median household\nincome ($)"))
        ) +
        theme(
          panel.grid.minor = element_blank(),
          panel.grid.major = element_blank(),
          strip.placement = "outside",
          axis.title.y = element_text(
            angle = 0,
            vjust = .5
          ),
          plot.margin = margin(7, 7, 7, 7),
          axis.line = element_line(),
          axis.line.y = element_line(),
          axis.ticks = element_line(),
          axis.text.y = element_text(vjust = .5, hjust = 1),
          plot.caption = element_text(
            size = rel(1),
            colour = "grey30"
          )
        ) +
        scale_y_continuous(
          labels = scales::percent_format(accuracy = 1),
          expand = expansion(mult = c(0, .05)),
          breaks = c(0, .15, .30, .45, .60)
        ) +
        scale_x_continuous(
          labels = scales::comma,
          expand = expansion(mult = c(0, .1))
        ) +
        labs(
          x = "", y = "Tree\ncanopy\n (%)",
          caption = # expression(italic(
          "Source: Analysis of Sentinel-2 satellite imagery (2021)\nand ACS 5-year estimates (2015-2019)" # ))
        )
fig_equity
ggsave("fig_equity.png",fig_equity,  width = 10, height = 5, units = "in", device = "png")
# ggsave("fig_equity.png",fig_equity,  width = 4, height = 5.5, units = "in", device = "png")

Hi income neighborhood:

hinc <- mn_bgs %>% select(GEOID) %>%
  rename(tract_string=GEOID) %>%
  right_join(df1 %>%
               filter(mdhhincnow > 100000,
                      canopy_percent > .35,
                      pop_density > 8,
                      pbipoc < .2))

pal <- colorNumeric(
  palette = "Greens",
  domain = hinc$canopy_percent)

leaflet(hinc) %>%
  addTiles() %>%
  addPolygons(color = "black",
              fillOpacity = .7,
              fillColor = ~pal(canopy_percent),
              popup = ~paste0(tract_string,
                              "<br>pop density: ", round(pop_density, 1), 
                              "<br>P bipoc: ", round(pbipoc(100), 2),
                              "<br>Tree canopy: ", round(canopy_percent, 2),
                              "<br>Median hhinc: $", prettyNum(mdhhincnow, big.mark = ","),
                              "<br>Temp: ", avg_temp,
                              "<br>NDVI: ", ndvi)) %>%
  addLegend("bottomright", 
            pal = pal,
            values = ~canopy_percent,
    opacity = 1
  ) 

Low income neighborhood

hinc <- mn_bgs %>% select(GEOID) %>%
  rename(tract_string=GEOID) %>%
  right_join(df1 %>%
               filter(mdhhincnow < 80000,
                      canopy_percent < .21,
                      pop_density > 8, 
                      pop_density < 15,
                      pbipoc > .4))

pal <- colorNumeric(
  palette = "Greens",
  domain = hinc$canopy_percent)

leaflet(hinc) %>%
  addTiles() %>%
  addPolygons(color = "black",
              fillOpacity = .7,
              fillColor = ~pal(canopy_percent),
              popup = ~paste0(tract_string,
                              "<br>pop density: ", round(pop_density, 1), 
                              "<br>P bipoc: ", round(pbipoc * 100, 2),
                              "<br>Tree canopy: ", round(canopy_percent, 2),
                              "<br>Median hhinc: $", prettyNum(mdhhincnow, big.mark = ","),
                              "<br>Temp: ", avg_temp,
                              "<br>NDVI: ", ndvi)) %>%
  addLegend("bottomright", 
            pal = pal,
            values = ~canopy_percent,
    opacity = 1
  ) 


# low = 270530001025 in camden
# high = 271230357002 in summit hill

bg_growingshade_data$ndvi[bg_growingshade_data$tract_string == "270530001025"] #96
bg_growingshade_data$avg_temp[bg_growingshade_data$tract_string == "271230357002"] #96

Temp

ndvilabs <- c(
      "<img src='./NDVI_.17.png' height='75' /><br>Low<br>green space",
      "<img src='./NDVI_.42.png' height='75' /><br>Moderate<br>green space",
      "<img src='./NDVI_.67.png' height='75' /><br>High<br>green space"
    )

df <- df1 %>%
        select(avg_temp, ndvi)

      tempplot <- ggplot(aes(x = ndvi, y = avg_temp), data = df) +
        geom_point(col = "grey40", alpha = .2, data = filter(df), na.rm = TRUE) +
        geom_smooth(method = "lm", formula = "y ~ x + I(x^2)", fill = NA, col = councilR::colors$councilBlue) +
        councilR::council_theme() +
        labs(
          x = "Amount of green space\n(maximum NDVI)", y = "Summer\nland surface\ntemperature\n(°F)",
          caption = "\nSource: Analysis of Sentinel-2 satellite imagery (2021)\nand Landsat 8 satellite imagery (2016)"
        ) +
        theme(
          panel.grid.minor = element_blank(),
          panel.grid.major = element_blank(),
          strip.placement = "outside",
          axis.title.y = element_text(
            angle = 0,
            vjust = .5
          ),
          plot.margin = margin(7, 7, 14, 7),
          axis.line = element_line(),
          axis.ticks = element_line(),
          axis.text.y = element_text(vjust = .5, hjust = 1),
          plot.caption = element_text(
            size = rel(1),
            colour = "grey30"
          ),
          axis.text.x.bottom = ggtext::element_markdown(size = 15)
        ) +
        scale_y_continuous(expand = expansion(mult = c(0, .05))) +
        scale_x_continuous(
          name = NULL,
          breaks = c(.17, .42, .67),
          labels = ndvilabs,
          position = "bottom"
        )
      tempplot
ggsave("fig_temp.png",tempplot,  width = 6, height = 5, units = "in", device = "png")

redline

rl <- bg_growingshade_data %>%
  select(tract_string, holc_pred, holc_pgrn, holc_pblu, holc_pylw, canopy_percent, ndvi, avg_temp) %>% #
  mutate_all(~replace(., . == 0, NA)) %>%
  mutate(flag = if_else(is.na(holc_pred) & is.na(holc_pblu) & is.na(holc_pgrn) & is.na(holc_pylw), "remove", "keep")) %>%
  filter(flag == "keep") %>% select(-flag) %>%
  pivot_longer(names_to = "holc", values_to = "percent", -c(tract_string, avg_temp, canopy_percent, ndvi)) %>%
  filter(!is.na(percent)) %>%
  mutate(grade = case_when(holc == "holc_pred" ~ "D (redlined)",
                           holc == "holc_pgrn" ~ "A (highest)",
                           holc == "holc_pylw" ~ "C",
                           holc == "holc_pblu" ~ "B"))
rl2 <- rl%>% filter(percent > .75) 

av <- rl2 %>% distinct(tract_string, .keep_all = TRUE) %>% summarise(avg_temp = mean(avg_temp))

redline_fig <- rl2 %>%
  ggplot(aes(y = fct_rev(grade), x = avg_temp, #(avg_temp - as.numeric(av)), 
             fill = (avg_temp))) +
  councilR::council_theme() +
  geom_vline(xintercept = as.numeric(av), color = "grey70") +
  # geom_vline(xintercept = 0, color = "grey70") +
  ggbeeswarm::geom_beeswarm(
            size = 2, 
            cex = 1.5,
            method = "compactswarm",
            na.rm = TRUE,
            pch = 21, color = "grey40"
          ) +
  scale_fill_distiller(palette = "RdBu") + 
  labs(x = "Land surface temperature during 2016 heat wave",
       #x = "2016 temperature difference from region average",
       y = "1934\nHome\nOwners'\nLoan\nCorporation\nrating",
       caption = "\nSource: Analysis of Landsat 8 satellite imagery (2016)\nand Equity Considerations dataset (2021) ") +
  theme(
          panel.grid.minor = element_blank(),
          panel.grid.major = element_blank(),
          strip.placement = "outside",
          axis.title.y = element_text(
            angle = 0,
            vjust = .5
          ),
          plot.margin = margin(7, 7, 7, 7),
          axis.line = element_line(),
          axis.line.y = element_line(),
          axis.ticks = element_line(),
          axis.text.y = element_text(vjust = .5, hjust = 1),
          plot.caption = element_text(
            size = rel(1),
            colour = "grey30"
          )
        ) +
  guides(fill = "none")
redline_fig
ggsave("fig_redline.png",redline_fig,  width = 7, height = 4.5, units = "in", device = "png")

Doc with canopy cover by community

library(tidyverse); library(sf)
ctu_list %>%
  st_drop_geometry() %>%
  dplyr::select(GEO_NAME, canopy_percent, min, max) %>%
  mutate(min = min/100,
         max = max / 100) %>%
  rename(`CTU Name` = "GEO_NAME",
         `Average tree cover (%)` = canopy_percent,
         `Lowest block group tree cover (%)` = min,
         `Highest block group tree cover (%)` = max) %>%
  write_csv("planit_trivia.csv")

carbon by tree stands from USForest Service EVALIDATOR

adapted from: https://extension.umn.edu/managing-woodlands/carbon-minnesota-trees-and-woodlands#manage-for-carbon-storage-2244060

access on web at: https://apps.fs.usda.gov/Evalidator/rest/Evalidator/fullreport?reptype=State&lat=0&lon=0&radius=0&snum=Aboveground and belowground carbon in live trees (at least 1 inch d.b.h./d.r.c), in short tons, on forest land&sdenom=Area of forest land, in acres&wc=272019&pselected=None&rselected=Forest Type MnDNR&cselected=Stand age 20 yr classes (0 to 100 plus)&ptime=Current&rtime=Current&ctime=Current&wf=&wnum=&wnumdenom=&FIAorRPA=FIADEF&outputFormat=HTML&estOnly=Y&schemaName=FS_FIADB.

usfs <- readxl::read_xlsx("../data-raw/evalidator fsusda.xlsx", 
                          skip = 1,
                          na = "-") %>%
  dplyr::select(-Total) %>%
  pivot_longer(names_to = "age",
               values_to = "carbon", 
               -`Forest Type MnDNR`) %>%
  rename("sps" = "Forest Type MnDNR") %>%
  filter(sps %not_in% c("Total", "Other",
                        "Other softwoods", "Non stocked",
                        "Cottonwood / Willow", "Eastern redcedar", "Balsam poplar", "White spruce", "Black spruce", "Balsam fir"))%>%
    mutate(age2 = case_when(age == "0-20 years" ~ "0-20",
                          age == "21-40 years" ~ "21-40",
                          age == "41-60 years" ~ "41-60",
                          age == "61-80 years" ~ "61-80",
                          age == "81-100 years" ~ "81-100",
                          age == "100+ years" ~ "100+",
                          TRUE ~ NA_character_)) %>%
  # group_by(sps) %>%
  mutate(age2 = factor(age2, levels=c("0-20",
                                    "21-40",
                                    "41-60",
                                    "61-80",
                                    "81-100",
                                    "100+")),
         label = case_when(age2 == "100+" ~ sps,
                           TRUE ~ NA_character_),
                  sps = fct_reorder(sps, desc(carbon))
         ) 


usfs_fig <- usfs %>%
  # group_by(sps)  %>%
  ggplot(aes(x = age2, y = carbon, color = sps, group = sps, shape = sps)) +
geom_line(lwd =1) +
  geom_point(aes(fill = sps),  size = 3) +
  # cowplot::theme_cowplot() +
  councilR::council_theme() +
  scale_color_brewer(palette = "Paired") +
  scale_fill_brewer(palette = "Paired") +
  # ggrepel::geom_label_repel(aes(label = label),
  #                 nudge_x = 1,
  #                 na.rm = TRUE
  #                 )
  ggrepel::geom_label_repel(
    aes("100+", carbon, label = label), col = "black", nudge_x = .5, direction = "y", hjust = "left",  ylim = c(-Inf, Inf),  xlim = c(NA, Inf)
  ) +
  scale_x_discrete(expand = expansion(mult = c(.05, 0.6))) +
  scale_shape_manual(values = rep((21:25), 3)) +
  guides(fill = "none", col = "none", shape = "none") +
  labs(x = "Stand age (years)", y = "Average\ncarbon\nstorage\n(US tons\nper acre)",
       caption = "US Forest Service EVALIDator v1.8.0.01") +
    theme(panel.grid.minor = element_blank(),
        panel.grid.major.x = element_blank(),
        axis.title.y = element_text(angle=0,
                                    vjust = .5),
        plot.margin = margin(7,7,7,7),
        legend.position = "none"#,
        # axis.ticks.x = element_line("grey")
        )

ggsave("./usfs_fig.jpg", usfs_fig,  width = 8, height = 5, units = "in", device = "jpg")

ampo

a <- bg_growingshade_main %>% 
  filter(variable %in% c("mdhhincnow", "canopy_percent", "pbipoc")) %>%
  dplyr::select(-name, -weights_scaled) %>%
  pivot_wider(names_from = variable, values_from = raw_value)

car::Anova(lm(canopy_percent ~ mdhhincnow, data = a))
car::Anova(lm(canopy_percent ~ pbipoc, data = a))

mn_bgs %>% sf::st_drop_geometry() %>% group_by(highest_priority) %>% count() %>%
  mutate(n = n/ nrow(mn_bgs))

Export data for storymaps

race_exp <- bg_growingshade_main %>%
  filter(variable == "pbipoc") %>%
  mutate(raw_value = raw_value * 100) %>%
  left_join(bg_geo %>% rename(bg_string = GEOID)) %>%
  dplyr::select(bg_string, raw_value, geometry) %>%
  rename(pbipoc = raw_value) %>%
  st_as_sf()

sf::st_write(race_exp, "/Users/escheh/Documents/GitHub/planting.shade/storymap-info/shapefiles/race.shp", append = FALSE)


canopy_exp <- bg_growingshade_main %>%
  filter(variable == "canopy_percent") %>%
  mutate(raw_value = raw_value * 100) %>%
  left_join(bg_geo %>% rename(bg_string = GEOID)) %>%
  dplyr::select(bg_string, raw_value, geometry) %>%
  rename(canopy_percent = raw_value) %>%
  st_as_sf()

sf::st_write(canopy_exp, "/Users/escheh/Documents/GitHub/planting.shade/storymap-info/shapefiles/canopy.shp", append = FALSE)


income_exp <- bg_growingshade_main %>%
  filter(variable == "mdhhincnow") %>%
  mutate(raw_value = raw_value * 100) %>%
  left_join(bg_geo %>% rename(bg_string = GEOID)) %>%
  dplyr::select(bg_string, raw_value, geometry) %>%
  rename(mdhhincnow = raw_value) %>%
  st_as_sf()

sf::st_write(income_exp, "/Users/escheh/Documents/GitHub/planting.shade/storymap-info/shapefiles/income.shp", append = FALSE)


copd_exp <- bg_growingshade_main %>%
  filter(variable == "COPD") %>%
  mutate(raw_value = raw_value * 100) %>%
  left_join(bg_geo %>% rename(bg_string = GEOID)) %>%
  dplyr::select(bg_string, raw_value, geometry) %>%
  rename(COPD = raw_value) %>%
  st_as_sf()

sf::st_write(copd_exp, "/Users/escheh/Documents/GitHub/planting.shade/storymap-info/shapefiles/copd.shp", append = FALSE)


highestp_exp <- bg_geo %>%
  right_join(mn_bgs %>%
               sf::st_drop_geometry()) %>%
               select(GEOID, highest_priority)
sf::st_write(highestp_exp, "/Users/escheh/Documents/GitHub/planting.shade/storymap-info/shapefiles/highestpriority.shp", append = FALSE)

# https://www.bls.gov/cpi/tables/supplemental-files/historical-cpi-u-202111.pdf
cpi05 <- 193.2
cpi21 <-  266.236
tribble(~`Benefits`, ~`Total ($)`, ~`SE ($)`, ~`$/tree`, ~`SE ($/tree)`, ~`$/capita`, ~`SE ($/capita)`,
 "Energy", 6824046, (483981), 34.36, (2.44), 8.79, (.62),
 "CO2", 826875, (58644), 4.16, (.3), 1.06, (.08),
 "Air quality", 1134334, (80450), 5.71, (.41), 1.46, (.1),
 "Stormwater", 9071809, (643399), 45.67, (3.24), 11.68, (.83),
 "Aesthetic/Other", 7076370, (501877), 35.63, (2.53), 9.11, (.65),
"Total Benefits", 24933434, (1766384), 125.53, (8.89), 32.10, (2.27)) %>%
  mutate(Total21 = `$/tree` * cpi21/cpi05) %>%
  ggplot(aes(x = Total21, y = Benefits)) +
  geom_point()

Export for PlanIt

ctu_list %>%
  sf::st_drop_geometry() %>%
  dplyr::select(GEO_NAME, canopy_percent, min, max) %>%
  mutate(min = min/100, 
         max = max/100) %>%
  rename("CTU Name" = GEO_NAME,
         "Average tree %" = canopy_percent,
         "Lowest block group tree %" = min,
         "Highest block group tree %" = max) %>%
  write_csv("./PlanIttriva.csv")


ctu_list %>%
  sf::st_drop_geometry() %>%
  dplyr::select(GEO_NAME, canopy_percent) %>%
  mutate(canopy_percent = case_when(canopy_percent < .10 ~ "<10%",
                                    canopy_percent <=.2 ~ "10-20%",
                                    canopy_percent <=.3~ "20-30%",
                                    canopy_percent <=.4 ~ "30-40%",
                                    canopy_percent <=.5 ~ "40-50%",
                                    TRUE ~ ">50%")) %>%
  count(canopy_percent) %>%
  mutate(total_com = sum(n),
         n = n / total_com) %>%
  mutate(canopy_percent = factor(canopy_percent, levels=c("<10%", "10-20%", "20-30%", "30-40%", "40-50%", ">50%"))) %>%
  ggplot(aes(y = fct_rev(canopy_percent), x = n,
             label = paste0(round(n*100,1), "%"))) +
  geom_bar(stat = "identity",
           width = .5)+
  labs(y = "2021\ntree\ncanopy", x = "Percent of communities") +
  theme_minimal() +
  theme(axis.title.y = element_text(angle = 0, vjust = .5)) +
  scale_x_continuous(labels = scales::percent, limits = c(0,.4)) +
  ggrepel::geom_label_repel(nudge_x = 10,segment.color = NA)

Metropolitan-Council/planting.shade documentation built on Feb. 25, 2024, 3:15 a.m.